  function [alph,bet,ct,xsim]=see_multbg(MU,SD,dens,mingr,stepgr,maxgr)
% function [alph,bet,ct,xsim]=see_multbg(MU,SD,dens,mingr,stepgr,maxgr)
% 
% Same as see_bg but can provide a vector of MU and SD, which will be
% plotted in the same graph. Will also print a comparison of percentiles, 
% automatically. 
% 
% Inputs
% ------
% MU    Mean  (must be < 1 for Beta)
% SD    Standard Deviation 
% dens  =1  or ='G' (string) for Gamma 
%       =2  or ='B' (string) for Beta 
% mingr  (optional)  Minimum value of the grid used to compute the pdf 
%                    (default is 0.01 for both densities)
% stepgr (optional)  Step of the grid used to compute the pdf  
%                    ( default is 0.01 )
% maxgr (optional)  Maximum value of the grid used to compute the pdf 
%                    (default is 1 for beta, 6 for Gamma)
%
% NOTES: 
%       1) If you want to change any of the three inputs of the gris you must 
%       call all three of them
%       2) Percentiles are not affected by these choices
%
% Output 
% [alph,beta]   location parameters  
% ct            cell array with the percentiles (also displayed in graph) 
% xsim          vector of 20,000 simulated values from a G(alph,bet) or 
%               Beta(alph,bet)  
%
% e.g.  see_bg(4,1,'G')  calls the function for the gamma with default
%       settings 
%
%       see_bg(0.5,0.2,1,0.07,0.001,0.8)  calls for a Beta and plots 
%       the interval 0.07 and 0.8
%
% alejandro justiniano 4/13/05 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Check all settings provided 
if MU < 0 | SD < 0 
    error('Mean and Variance must be positive for these densities'); 
end

if nargin < 4 
    disp('Default grid settings will be used') 
    mingr=[]; maxgr=[]; stepgr=[]; 
elseif nargin >= 4 & nargin < 6 
    error('Need to provide Min,Step, Max of the grid as inputs') 
end 

if isstr( dens ) == 1 
    dens=upper(dens);
    if dens =='G'
        dens=1 
    elseif dens=='B'
        dens=2; 
    else 
        error('Density must be G or 1, B or 2') 
    end
else 
    if dens < 1 | dens > 2 
        error('Density must be G or 1, B or 2')
    end 
end 

switch dens 
    case 1 
        bet=(SD^2)/MU; 
        alph=MU/bet; 
        [ct,xsim]=plot_gamma(alph,bet,mingr,stepgr,maxgr);
    case 2
        if MU > 1; error('Mean must be less than 1'); end;
        bet=MU*( (1-MU)^2 )/(SD^2) - ( 1 - MU );
        alph=MU*bet/(1-MU); 
        if bet < 0 | alph < 0; 
            error('Check the settings of your beta') 
        end 
        [ct,xsim]=plot_beta(alph,bet,mingr,stepgr,maxgr); 
end    
